% fig3b.m
clear; clc; close all;

%% 1. 参数设置
mu = -0.1;
x0 = 0.8;     
q0 = 0.6;
k_start = 1.25;
k_end   = 2.0;
dk      = 0.002;
k_vals  = k_start : dk : k_end;

% 迭代相关设置
transient_steps = 300;  % 舍去初始瞬态步数
record_steps    = 100;  % 用于分岔图记录的步数

% 用于存放计算结果
LE1_array = zeros(size(k_vals));
LE2_array = zeros(size(k_vals));

x_bif = [];  % 存储分岔图 x
q_bif = [];  % 存储分岔图 q
k_bif = [];  % 对应的 k

index = 1;
for k = k_vals
    % 每次从同一初值开始
    x = x0;
    q = q0;
    
    % (a) 排除瞬态
    for i = 1:transient_steps
        [x, q] = MSM(x, q, mu, k);
    end
    
    % (b) 记录 record_steps 步作为分岔数据
    x_record = zeros(1, record_steps);
    q_record = zeros(1, record_steps);
    x_record(1) = x;
    q_record(1) = q;
    for i = 2:record_steps
        [x_record(i), q_record(i)] = MSM(x_record(i-1), q_record(i-1), mu, k);
    end
    
    x_bif = [x_bif, x_record];
    q_bif = [q_bif, q_record];
    k_bif = [k_bif, k * ones(1, record_steps)];
    
    % (c) 计算有限时间李雅普诺夫指数
    [LE1, LE2] = LEs_MSM(mu, k, x0, q0);
    LE1_array(index) = LE1;
    LE2_array(index) = LE2;
    index = index + 1;
end

%% 绘图
figure('Name','Fig3(b) - 2D-MSM');
% 上子图：LEs 曲线
subplot(2,1,1);
hold on; grid on;
plot(k_vals, LE1_array, 'r', 'LineWidth',1.5, 'DisplayName','LE_1');
plot(k_vals, LE2_array, 'b', 'LineWidth',1.5, 'DisplayName','LE_2');
yline(0,'k--','LineWidth',1);  
xlabel('k'); ylabel('LEs');
title(['2D-MSM: \mu=', num2str(mu), ', (x_0,q_0)=(', num2str(x0), ',', num2str(q0),')']);
legend('Location','Best');

% 下子图：分岔图
subplot(2,1,2);
hold on; grid on;
scatter(k_bif, x_bif, 1, 'g', 'filled', 'DisplayName','x');
scatter(k_bif, q_bif, 1, 'm', 'filled', 'DisplayName','q');
xlabel('k'); ylabel('x, q');
legend('Location','Best');

%% 子函数
function [x_next, q_next] = MSM(x, q, mu, k)
    % 2D-MSM 映射:
    % x_{n+1} = mu*sin(2*pi*x_n) + k*cos(q_n)*x_n
    % q_{n+1} = q_n + x_n
    x_next = mu*sin(2*pi*x) + k*cos(q)*x;
    q_next = q + x;
end

function [LE1, LE2] = LEs_MSM(mu, k, x0, q0)
    % 采用 QR 正交化方法计算有限时间李雅普诺夫指数
    MaxIter = 1000;
    Transient = 200;
    
    x = x0; q = q0;
    e1 = [1; 0];
    e2 = [0; 1];
    sum_log1 = 0; sum_log2 = 0;
    
    for i = 1:MaxIter
        % 雅可比矩阵计算
        J = [ mu*(2*pi)*cos(2*pi*x) + k*cos(q), -k*sin(q)*x;
              1,                             1 ];
        % QR 分解
        e1_new = J * e1;
        e2_new = J * e2;
        [Q, R] = qr([e1_new e2_new], 0);
        e1 = Q(:,1);
        e2 = Q(:,2);
        
        if i > Transient
            sum_log1 = sum_log1 + log(abs(R(1,1)));
            sum_log2 = sum_log2 + log(abs(R(2,2)));
        end
        
        [x, q] = MSM(x, q, mu, k);
    end
    
    eff_iter = MaxIter - Transient;
    LE1 = sum_log1 / eff_iter;
    LE2 = sum_log2 / eff_iter;
end
